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Abstract 

Background: Bithynia siamensis goniomphalos is the snail intermediate host of the liver fluke, Opisthorchis viverrini, the 
leading cause of cholangiocarcinoma (CCA) in the Greater Mekong sub-region of Thailand. Despite the severe public health 
impact of Op/sf/?o/rh/s-induced CCA, knowledge of the molecular interactions occurring between the parasite and its snail 
intermediate host is scant. The examination of differences in gene expression profiling between uninfected and 0. viverrini- 
infected B. siamensis goniomphalos could provide clues on fundamental pathways involved in the regulation of snail- 
parasite interplay. 

Methodology/Principal Findings: Using high-throughput (lllumina) sequencing and extensive bioinformatic analyses, we 
characterized the transcriptomes of uninfected and 0. viverrin /-infected B. siamensis goniomphalos. Comparative analyses of 
gene expression profiling allowed the identification of 7,655 differentially expressed genes (DEGs), associated to 43 distinct 
biological pathways, including pathways associated with immune defense mechanisms against parasites. Amongst the 
DEGs with immune functions, transcripts encoding distinct proteases displayed the highest down-regulation in Bithynia 
specimens infected by 0. viverrini; conversely, transcription of genes encoding heat-shock proteins and actins was 
significantly up-regulated in parasite-infected snails when compared to the uninfected counterparts. 

Conclusions/Significance: The present study lays the foundation for functional studies of genes and gene products 
potentially involved in immune-molecular mechanisms implicated in the ability of the parasite to successfully colonize its 
snail intermediate host. The annotated dataset provided herein represents a ready-to-use molecular resource for the 
discovery of molecular pathways underlying susceptibility and resistance mechanisms of B. siamensis goniomphalos to 0. 
viverrini and for comparative analyses with pulmonate snail intermediate hosts of other platyhelminths including 
schistosomes. 
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Introduction 

Cholangiocarcinoma (CCA, or bile-duct cancer) is one of the 
most common malignant cancers in Southeast Asia, with incidence 
rates ranging from 93.8 to 317.6 per 100,000 people/year in some 
districts of Northeast Thailand alone [1,2]. Chronic infections by 
the food-borne human liver fluke Opisthorchis viverrini have been 
associated with the occurrence of CCA in several Southeast Asian 
countries (including Thailand and Laos) [2,3], which has led to the 
inclusion of this parasite amongst the Group 1 carcinogens by the 



International Agency for Research on Cancer [4]. Typically, CCA 
caused by 0. viverrini is diagnosed 30-40 years after primary 
infection, and death occurs within 3-6 months post-diagnosis [5] . 
The life cycle of 0. viverrini is complex, involving two intermediate 
hosts and multiple developmental stages [6]. Briefly, piscivorous 
definitive hosts, including humans, shed the embryonated eggs 
through the faeces; in aquatic environments, eggs are ingested by 
freshwater prosobranch snails of the genus Bithynia, in which the 
motile embryos ( = miracidia) hatch and develop into sporocysts. 
In the body of the first intermediate hosts, the sporocysts undergo 
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Author Summary 

Despite recent significant advances in knowledge of the 
fundamental biology of the carcinogenic liver fluke 
Opisthorchis viverrini, little is known of the complement 
of molecular interactions occurring between this parasite 
and its prosobranch snail intermediate host, Bithynia 
siamensis goniomphalos. The determination of such inter- 
actions is a key, necessary component of the development 
of future integrated control strategies for liver fluke 
infection and associated bile duct cancer. Here, we use 
cutting-edge high-throughput sequencing technologies 
and advanced bioinformatic analyses to characterize, for 
the first time, qualitative and quantitative differences in 
gene expression between uninfected and O. viverrini- 
infected B. siamensis goniomphalos collected from an 
endemic region of Northeast Thailand. The analyses led to 
the identification of a number of molecules putatively 
involved in immune defense pathways against invading 0. 
viverrini, and of key biological mechanisms potentially 
implicated in the ability of the parasite to successfully 
colonize its snail intermediate host. We believe that this 
ready-to-use molecular resource will provide the scientific 
community with new tools for the development of 
strategies to control the spread of liver fluke infection 
and the resulting bile duct cancer. 

asexual reproduction through the stages of rediae and cercariae. 
Subsequently, the motile cercariae are released into the aquatic 
environment where they actively seek the second intermediate 
hosts, i.e. cyprinid fishes (e.g. Puntius spp.). Upon penetrating the 
skin of these fishes, the parasites encyst as metacercariae under the 
skin and/or in muscular tissue. Humans acquire the infection by 
ingesting raw or undercooked infected fishes [5,6]. Following 
ingestion, the metacercariae excyst in the duodenum and migrate 
as juvenile flukes to the intra-hepatic bile ducts, where they 
develop to adult hermaphrodite flukes within ~4 weeks. The 
sexually mature adults release eggs that, via the bile, are excreted 
with the faeces into the environment [6]. 

Over the last decade, much effort has been directed towards the 
elucidation of the molecular mechanisms that determine the 
occurrence of CCA in 0. viverrini-mitcted individuals [2,7-9], and 
advances in high-throughput sequencing technologies and bioin- 
formatics are contributing towards a better understanding of the 
systems biology of this parasite [10,1 1]. While these advancements 
undoubtedly provide a foundation for the development of novel 
strategies for the diagnosis, treatment and prevention of 0. viverrini- 
induced CCA [3], long-term control of this devastating disease 
stricdy relies on the development of integrated approaches, 
targeting the parasite as well as its intermediate hosts [12,13]. 
Thus far, large-scale molecular investigations of snail species 
acting as intermediate hosts of human parasites have mainly 
involved pulmonate species of the genus Biomphalaria, which 
harbor Schistosoma blood flukes (see Biomphalaria glabrata Genome 
Initiative at http:/ /biology.unm.edu/biomphalaria-genome/index. 
html; [14,15]) while, until recentiy, knowledge of the fundamental 
molecular biology of prosobranch snails has remained limited 
[12,13]. Our group has recently utilized high-throughput (Illumina) 
sequencing and advanced bioinformatic analyses of sequence data 
to generate and characterize the first reference transcriptome of 
Bithynia siamensis goniomphalos (Gastropoda, Bithyniidae), the snail 
intermediate host of 0. viverrini in areas of Northeast Thailand, Lao 
PDR, Cambodia and South Vietnam [13]. This sequence dataset 
represents a valid resource for molecular studies of these organisms 
and provides a solid basis for investigations of snail-parasite 



interactions at the molecular level. In the present study, we employ 
high-throughput (Illumina) sequencing and bioinformatic technol- 
ogies to characterize differences in levels of gene transcription 
between uninfected and O. viverrini-iaStcted B. siamensis goniomphalos 
and identify key molecules and biological pathways inferred to be 
associated with the processes of invasion and survival of the parasite 
in the snail intermediate host. 

Materials and Methods 

Procurement of snail materials and assessment of 
infection status 

Adult B. siamensis goniomphalos snails were collected from natural 
freshwater bodies in Amphoe Muang, Khon Kaen Province, 
Thailand. The soft tissues of each snail were removed from the 
shells and washed three times with sterile distilled water. A portion 
of the hepatopancreas was isolated from each snail and subjected 
to DNA extraction using the Direct PCR kit (Thermo Scientific) 
according to manufacturer's instructions, while the remaining 
material was preserved in RNAlater (Sigma) for subsequent 
extraction of total RNA. A PCR targeting the 0. viverrini pOV-A6 
gene was employed to screen individual snails for infection by the 
liver fluke following the amplification protocol described by 
Sermswan et al. [16]. Amplicons were visualized on a 1.5% 
agarose gel stained in ethidium bromide and the presence of bands 
of approximately 330 bp in size was considered diagnostic for 0. 
viverrini infection [17]. 

RNA isolation and Illumina sequencing 

Two distinct pools of soft tissues from uninfected (UB; n = 4) 
and O. viverrini-'mtected snails (IB; n = 4) were subjected to isolation 
of total RNA using the Trizol reagent (Invitrogen), followed by 
DNAse I treatment (Promega) according to the manufacturer's 
instructions. The amounts and integrity of total RNA were 
determined using a 2100 BioAnalyzer (Agilent Technologies). 
Polyadenylated (PolyA+) RNA was purified from 10 (Xg of total 
RNA from each pool of uninfected and 0. viverrini-infected snails 
using Sera-Mag oligo(dT) beads, fragmented to a length of 1 00— 
500 nucleotides and reverse transcribed to cDNA using random 
hexamers. The size-fractionated cDNA was end-repaired and 
adaptor-ligated according to the manufacturer's protocol (Illu- 
mina). Ligated products of ~300 bp were excised from agarose 
gels and PCR-amplified (15 cycles) [18]. Products were cleaned 
using a MinElute column (Qiagen, Netherlands) and paired-end 
sequenced on an Blumina HiSeq 2000 [19] according to the 
manufacturer's protocol. 

Bioinformatic analyses 

Following removal of adapter sequences and sequences with 
suboptimal read quality (i.e., PHRED score of 32.0), the 
remaining 90 bp single-read sequences generated from the non- 
normalized cDNA libraries from UB and IB, as well as the raw 
sequence data generated previously [13] were pooled and 
assembled using the program Trinity [20] (http://trinityrnaseq. 
sourceforge.net/). Briefly, a representative set of transcripts was 
assembled using the software Inchworm; next, the Chrysalis 
software was used to cluster portions of alternatively spliced 
transcripts and/or unique portions of paralogous genes and 
construct a de Brujin graph for each cluster of transcripts. Finally, 
each de Brujin graph was analysed and alternatively spliced 
isoforms and transcripts derived from paralogous genes were 
resolved using the Butterfly software [20] (http://trinityrnaseq. 
sourceforge.net/). The non-redundant, assembled dataset was 
subsequently compared to transcriptomic sequence data available 
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for 0. viverrini [10,11] as well as the entire genome sequences of 
Clonorchis sinensis, S. mansoni, S.japonicum and S. haematobium [10,21- 
24]; sequences with high similarity to known 0. viverrini and/or 
other molecules from parasitic trematodes (e-value cut-off: le-15) 
were excluded [25]. The non-redundant transcriptomic dataset for 
B. siamensis goniomphalos was then annotated using an established 
approach [23,25]. Briefly, assembled contigs were compared using 
the BLASTn and BLASTx algorithms [26] to sequences available 
in the nucleotide sequence collection (nt) of NCBI (www.ncbi.nlm. 
nih.gov), and in the non-redundant (nr) (www.ncbi.nlm.nih.gov), 
Swiss-Prot (http://expasy.org/), Kyoto Encyclopedia of Genes 
and Genomes (KEGG; http:/ /www.genome.jp/kegg/) and Clus- 
ters of Orthologous Groups (COG; http://www.ncbi.nlm.nih. 
gov/COG/) databases, respectively (e-value<0. 00001), to identify 
putative homologues in other mollusks and organisms other than 
mollusks. In addition, the non-redundant B. siamensis goniomphalos 
transcriptome was compared, using a BLASTn algorithm, to an 
in-house built nucleotide sequence database containing 86,837 
expressed sequence tags (ESTs) available for B. glabrata (http:// 
www.ncbi.nlm.nih.gov/ nucest/?term = biomphalaria+glabrata). 

Proteins with the highest sequence similarity to conceptually 
translated B. siamensis goniomphalos transcripts were assigned Gene 
Ontology terms (GO; http://www.geneontology.org/) [27], ac- 
cording to the categories 'Biological Process', 'Cellular Compo- 
nent' and 'Molecular Function', using the Blast2GO software [28]. 
B. siamensis goniomphalos contigs with no homologues in available 
protein databases were conceptually translated using the ESTScan 
software and annotated using the same approach [29]. 

The short reads from each UB and IB dataset were then 
separately mapped to the non-redundant sequence data using 
SOAPAligner, such that each read was mapped to a unique 
transcript with a minimum alignment overlap of 40 nucleotides 
and a maximum of three mismatches per read. Reads that mapped 
to more than one transcript ("multi-reads") were randomly 
allocated to a unique transcript, such that they were mapped 
only once. Comparative analyses of gene transcription levels 
between UB and IB were performed using the fragments per 
kilobase of exon per million fragments method (i.e., FPKM) [30]. 
The analysis of statistical difference of transcription was deter- 
mined using a method developed for serial analysis of gene 
expression (SAGE) and applied to RNA-seq data [31]. Statistical 
significance was set at & p value of £0.01 and, to control for errors 
associated with multiple pairwise comparisons, a false-discovery 
rate (FDR) correction was applied to the data set [32]. GO and 
KEGG functional enrichment analyses for B. siamensis goniomphalos 
genes whose transcription levels were altered by infection with the 
parasite, were performed using the software GO::TermFinder with 
multiple testing correction using the Benjamini & Hochberg 
calculations of FDR (cut-off: FDR<0.01) [33,34]. 

Results 

A total number of 60,210,044 and 56,426,082 raw reads were 
generated from UB and IB, respectively; of these, 55,137,440 
(92%) and 52,476,468 (92%), respectively, were retained for 
further analyses (Table 1). All single reads generated in the present 
study have been deposited in the Sequence Read Archive (SRA) 
database of NCBI (http://www.ncbi.nlm.nih.gov/sra) under 
accession numbers SRR1046838 (UB) and SRR1046839 (IB), 
respectively. The assembly yielded 43,259 contigs, of which 20,171 
displayed significant sequence similarity (e-value cut-off: le-15) to 
publicly available sequences from 0. viverrini, C. sinensis, S. mansoni, 
S. japonicum and/or S. haematobium (data not shown) and were 
therefore excluded from subsequent analyses. All remaining 



23,088 contigs (average length = 775 nt±679; maximum 
length = 12,119 nt; minimum length = 200; N50 = 671; cf. 
Table 1), matched previously assembled sequence data for B. 
siamensis goniomphalos, and displayed significant sequence similarity 
(e-value cut-off: le-05) with available ESTs from B. glabrata (not 
shown). Of the assembled sequence dataset, 5,461 (23.6%) and 
13,283 (57.5%) contigs had homologues in the non-redundant 
nucleotide (nt) and protein (nr) sequence databases of NCBI, 
respectively (e-value cut-off: le-05; Table 1). Of the conceptually 
translated contigs with significant similarity to sequences in the nr 
database, 18.4% matched known proteins from the bivalve 
mollusk Crassostrea gigas (Table SI). This Transcriptome Shotgun 
Assembly Project has been deposited at DDJB/EMBL/ GenBank 
under the accession GAQOJ30000000. 

Of the 17,849 contigs that encoded a predicted protein (cf. 
Table 1), 10,761 (60.3%) had homologues in the Swiss-Prot 
database (cf. Table 1); these included mucins (n = 378), histones 
(n=318), collagens (n = 251) and serine/threonine protein phos- 
phatases (STPs; n = 181) and kinases (STKs; n = 119) (Table S 1). A 
total of 5,950 predicted proteins (33.3%) could be assigned GO 
terms, while 9,497 (53.2%) matched homologous proteins in the 
KEGG database associated to 259 distinct biological pathways 
(Table 1). The 5,287 (29.6%) predicted proteins with matches in 
the COG database could be assigned to at least one of 25 
functional categories, of which 'general function prediction' (84%), 
'translation, ribosomal structure and biogenesis' (50.2%) and 'post- 
translational modification, protein turnover and chaperones' 
(41%) were the most represented (Tables 1 and SI). 

Comparative analyses of levels of gene expression revealed 
6,920 and 735 DEGs in UB and IB, respectively (FPKM>2.0; 
Table 1, Figure 1). Of the entire complement of molecules up- 
regulated in UB, 4,488 (64%) did not match known proteins in 
the nr database (Table S2); of the remaining 2,462, 7% were 
mucins, followed by uncharacterized proteins (2%) and immu- 
noglobulin receptor homologues (0.6%). Of the 208 DEGs in IB 
with known homologues in the nr database (cf. Table S2), 38% 
were uncharacterized proteins, followed by actins (1%) and 
transcription factors (0.7%) (Table S2). Two haemocyanins and a 
Kunitz-type protease inhibitor were amongst the most signifi- 
cantly up-regulated DEGs in UB, whereas a histone and dynein 
displayed the highest expression levels in IB (Table 2). GO- 
enrichment analysis of DEGs in UB and IB revealed significant 
enrichment of the GO terms 'cellular process' (80.8% of DEGs), 
'organic substance metabolic process' (53.3%) and 'cellular 
metabolic process' (52.1%) within 'biological process', 'cell' 
(90.4%), 'intracellular part' (85.4%) and 'organelle' (68.7%) 
within 'cellular component' and 'protein binding' (28.3%), 
'hydrolase activity' (25.1%) and 'ribonucleoside binding' within 
'molecular function' (Figure 2a). KEGG pathway enrichment 
analysis identified a total of 43 pathways linked to DEGs, 
including 'RNA transport' (8.6%), 'Huntington's disease' (7.1%) 
and 'spliceosome' (8.5%) (Figure 2b). Out of 375 DEGs predicted 
to be involved in biological pathways linked to immune responses 
(cf. Figure 2b), homologues of metalloproteases, and actins and 
heat shock proteins displayed the largest up-regulation in UB and 
IB, respectively (Figure 3). 

Discussion 

Despite the high prevalence of infection with 0. viverrini in 
humans in the Greater Mekong sub-region (i.e. >90% in lowland 
districts of Lao PDR), and the large number of cyprinid fishes 
which may act as potential second intermediate hosts for this 
parasite [35-37], snail infection rates are extremely low, i.e. 
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Table 1. Summary of the RNA-Seq data from uninfected and Opisthorchis v/Verrin/'-infected Bithynia siamensis goniomphalos snails 
(UB and IB, respectively), prior and following assembly, and results of bioinformatic analyses and annotation. 




Raw reads (paired-end) 


60,210,044 (UB); 56,426,082 (IB) 


Assembly 


Contigs (average length ± SD; min-max length) 


23,088(775±679) 


N50 


671 


GC content (%) 


46.32 


Annotation 


Contigs containing an Open Reading frame (%) 


1 7,849 (77.3) 


Inferred via BLAST alignment 


16,398 


ESTScan predictions 


1,451 


Returning SwissProt hits (%) 


10,761 (60.3) 


Gene Ontology (%) 


5,950 (33.3) 


Number of parental Biological Process terms 


23 


Cellular component 


16 


Molecular function 


17 


Returning a KEGG result (%) 


9,497 (53.2) 


Number of predicted biological pathways 


259 


Returning a COG result (%) 


5,287 (29.6) 


Number of predicted functional categories 


25 


DEGs 


Total number of DEGs (%) 


7,655 (33) 


Up-regulated in UB 


735 


Up-regulated in IB 


6,920 


The number of statistically significant (p<0.01) differentially expressed genes (DEGs) in each UB and IB is also reported. 
doi:1 0.1 371 /journal.pntd.0002765.t001 



typically <1% [37-41] with peaks of ~ 3 % recorded in wetlands of 
Thailand and Lao PDR [42]. These observations lead to the 
hypothesis that infections of Bithynia snails by O. viverrini may 
stimulate a cascade of immune-molecular events similar to those 
induced in Biomphalaria spp. in response to the invasion by 
schistosomes, which are aimed at eliminating and/or preventing 
the dissemination of the infection in the snail intermediate hosts 
[14,15,43-45]. While knowledge of the complex mechanisms 
underlying the activation of immune defense pathways in 
Biomphalaria snails is expanding [46,47], studies of the immune- 
molecular interactions in the Bithynia- Opisthorchis system are 
lacking. By undertaking the first large-scale comparative analysis 
of gene expression in UB and IB, we hoped to provide the 
community with a ready-to-use molecular infrastructure for in- 
depth studies of biological pathways specifically involved in snail- 
parasite interactions. An obvious limitation of this study was 
represented by the fact that the snail specimens subjected to RNA- 
Seq were collected from natural freshwater bodies, and it was 
therefore impossible to establish whether the uninfected Bithynia 
snails had been previously exposed to infection by 0. viverrini and/ 
or other parasites, or whether a previous infection had successfully 
been cleared. However, unlike the monoecious Biomphalaria spp., 
the dioecious nature of Bithynia poses some challenges for the 
rearing and maintenance of this snail species under experimental 
conditions. Indeed, it has been reported that Bithynia spp. 
specimens reared in the laboratory are smaller and morpholog- 
ically altered when compared to the wild counterparts [40,48]; 
since changes in body morphology are likely to correspond 
to alterations in gene regulation and expression, both in the 
steady-state and following infection by O. viverrini, we preferred to 



investigate transcriptional differences between UB and IB 
collected in the field. 

Despite displaying significant sequence similarities to ESTs 
from B. glahrata, a significant proportion of B. siamensis goniomphalos 
DEGs (65%) could not be functionally annotated based on 
homology to known protein sequence data available in public 
databases. This finding is likely to be related to the limited 
number of annotated proteins available from mollusks [49], 
which impairs meaningful functional analysis of large-scale 
transcriptomic datasets. For DEGs matching available protein 
sequences, GO and KEGG enrichment analyses identified a 
range of biological pathways linked to immune responses, e.g. T 
and B cell signaling pathways (cf. Figure 2b), consistent with the 
results of previous investigations of levels of gene transcription in 
Biomphalaria snails exposed to infection by schistosomes [14,43- 
45,50,51]. Amongst the immunity-associated DEGs, two prote- 
ases, a metalloprotease and a cysteine protease (cathepsin L) 
displayed the largest degree of up-regulation in UB when 
compared with IB (cf. Figure 3). The role of proteolytic enzymes, 
such as cathepsins, in the defense of mollusks against invading 
trematodes (e.g. in the breakdown of components of invading 
parasites) has been well documented [52,53]. For instance, 
enhanced proteolytic activity has been described from extracts 
from S. maw-sora-resistant compared to susceptible B. glahrata [53], 
while other studies have consistendy detected transcripts encod- 
ing proteases amongst the molecules up-regulated in schistosome- 
resistant snail lines [14,44]. In particular, two distinct expressed 
sequence tags (ESTs) encoding cathepsins L had been previously 
identified as highly expressed in B. glahrata snails susceptible to 
infection by the trematode parasite Echinostoma caproni [54] and 
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■ Up-regulated 

■ Down-regulated 




-4 -3 -2 -1 0 1 2 3 4 5 6 

loglO(FPKM) - IB 



Figure 1 . Scatter plot comparing log2 ratios of FPKM expression values for uninfected and Opisthorchis wVemm-infected Bithynia 
siamensis goniomphalos snails (UB and IB, respectively). 

doi:1 0.1 371 /journal.pntd.0002765.g001 



resistant to S. mansoni [14], respectively. Based on the results of 
bioinformatic analyses of these two ESTs, Lockyer et al. [14] 
hypothesized that these two sequences might have represented 
two non-overlapping fragments of the same gene, and that the 
discrepancies in gene expression which had been observed 
followed infection by E. caproni and S. mansoni might have been 
due to "the different response elicited by two different parasite 
species" [14]. Comparative analyses between the abovemen- 
tioned Biomphalaria ESTs and the cathepsin L-encoding sequence 
identified in the present study suggest that the latter transcript is 
indeed distinct from the two previously identified (data available 
from the corresponding author upon request). Moreover, in 
addition to the transcript encoding a 'novel' cathepsin L, a 
molecule encoding a ficolin was identified amongst the most 
highly expressed transcripts in UB (cf. Figure 3). Ficolins are a 
group of lectin-like defense proteins that, in vertebrates, are 
known to recognize carbohydrate molecules on pathogens and/ or 
apoptotic and necrotic cells [55]. Structurally, ficolins consist of a 
small N-terminal domain, a collagen-like domain and a 
fibrinogen-like domain similar to the C-terminal domains of the 
beta and gamma chains of fibrinogen [56]. These molecules, 
which belong to a large protein family known as Fibrinogen- 
Related Proteins (FREPs), play key roles in the innate immunity 
of both vertebrate and invertebrate organisms [5 7] . In particular, 
in mollusks, FREPs are believed to exert 'antibody-like' 
properties in biological mechanisms of recognition and binding 



of invading pathogens [58,59]. In B. glabrata, transcription of 
FREP-encoding genes is enhanced after exposure of the snail to 
both S. mansoni and E. paraensei [43,60,61]; furthermore, FREPs 
have been shown to bind sporocysts of both S. mansoni and E. 
paraensei, thus actively participating in molecular events linked to 
defense against pathogens [62] . A mucin and an immunoglobulin 
receptor homologue were also amongst highly transcribed 
molecules in UB (cf. Figure 3). While the up-regulation of 
transcripts encoding mucins could be attributed to a role of these 
proteins in the establishment of a physical defense barrier against 
invading pathogens [63], the detection of an immune receptor- 
like molecule amongst the DEGs is interesting, since mollusks 
lack adaptive immunity. However, in a previous study, Bayne and 
co-workers [64] demonstrated that hemocytes of B. glabrata were 
unable to recognize IgG-bound antigens of S. mansoni sporocysts, 
thus providing evidence for the existence of specific antigen- 
binding sites on snail hemocytes [64]. Results from the present 
study support this hypothesis. 

While it is tempting to speculate that the identification of 
homologues of a protease, a FREP and an immunoglobulin- 
receptor amongst the most highly expressed transcripts in UB 
(collected from the same natural water bodies as IB) may be 
associated with a 'successful' defense response against 0. viverrini, 
which ultimately resulted in prevention or clearance of the 
infection (similarly to resistant strains of B. glabrata), such a 
hypothesis requires rigorous testing. 
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Table 2. The top-ten up-regulated genes in uninfected and Opisthorchis viverrin /-infected Bithynia siamensis goniomphalos snails 
(UB and IB, respectively). 





Contig ID 


Size (bp) 


FPKM 


Top match in nr* database 


ko definition** 


Top match in Biomphalaria 
glabrata*** 


UB 


CL4735.Contig4 


710 


19.50 


Haemocyanin [Aplypsio californica] 


Haemocyanin 


Tyrosinase-domain containing 
protein 


CL3120.Contig5 


257 


19.08 


Uncharacterised protein 


Uncharacterised 


Fibrinogen-related domain 
containing protein 


CL1900.Contig1 


534 


17.74 


Keyhole limpet haemocyanin 
[Megathura crenulata] 


Haemocyanin 


Tyrosinase-domain containing 
protein 


CL303.Contig1 


683 


17.54 


Uncharacterised protein 


Uncharacterised 


Hypothetical protein 


CL8571.Contig1 


417 


17.47 


Kyelin/chordin-like protein 
[Crassostrea gigas] 


Kunitz-type protease inhibitor 


Kunitz-type protease inhibitor 


CL1713.Contig1 


1159 


17.26 


Hypothetical protein 
[Branchiostoma fioridae] 


Lectin 


C-type lectin domain containing 
protein 


CL6142.Contig2 


358 


17.21 


Uncharacterised protein 


Uncharacterised 


Hypothetical protein 


CL10409.Contig1 


328 


17.17 


Uncharacterised protein 


Uncharacterised 


Hypothetical protein 


CL800.Contig6 


1116 


17.13 


Hypothetical protein 
[Branchiostoma fioridae] 


Fibropellin 


Chitin-binding, multiple EGF and 
TSP domain-containing protein 


CL8747.Contig1 


1315 


17.10 


Hypothetical protein 
[Capitella teleta] 


Fructose-biphosphatase 
aldolase 


Fructose-biphosphate aldolase 


IB 


CL10219.Contig4 


583 


17.55 


Histone [Harpegnathos saltator] 


Histone 


Histone 


CL636.Contig4 


287 


16.44 


Uncharacterised protein 


Uncharacterised 


Hypothetical protein 


CL8536.Contig2 


282 


16.41 


Uncharacterised protein 


Uncharacterised 


Hypothetical protein 


CL1651.Contig1 


388 


16.21 


Wiskott-Aldrich syndrome 
protein [Myotis davidii] 


Splicing factor 


Thioester-containing protein 


CL4675.Contig3 


930 


15.97 


ADP-ribosylation factor 
[Wuchereria bancroftf] 


ADP-ribosylation factor 


ADP-ribosylation factor 


CL4639.Contig1 


474 


16.10 


Diaphanous-like protein 
[Pteropus alecto] 


Formin 


Relish protein 


CL8803.Contig4 


1103 


16.07 


Hypothetical protein 
[Branchiostoma fioridae] 


Gamma-aminobutyric 

acid receptor-associated protein 


Hypothetical protein 


CL9377.Contig2 


726 


16.06 


Uncharacterised protein 


Uncharacterised 


Calcium-binding protein 


CL1031.Contig2 


753 


15.98 


Uncharacterised protein 


Uncharacterised 


p53-like transcription factor 


CL43.Contig1 


345 


15.76 


Dynein light-chain 
[Crassostrea gigas] 


Dynein 


Dynein 



*non-redundant database. 

**KEGG-Orthology based annotation definition. 

***Based on comparisons with the preliminary Biomphalaria glabrata genome assembly v4.3 (http://biology.unm.edu/biomphalaria-genome/). 
doi:1 0.1 371/joumal.pntd.0002765.t002 



Amongst the DEGs that could be linked to pathways associated 
with immune responses, actins and heat-shock proteins (HSPs) 
were identified amongst the molecules with the largest degree of 
up-regulation in IB (cf. Figure 3). The function/ s of HSPs in stress 
responses of hemocytes of B. glabrata in response to schistosome 
infection has been documented [65,66]; conversely; the role/s 
played by actins in snail-parasite relationships is unclear. Since 
defense against pathogens in mollusks mainly relies on the 
migration of circulating hemocytes to the sites of infection and 
tissue damage, it could be hypothesized that the relative increase 
in transcription of actin-encoding genes might reflect a general 
increase in cellular motility, as well as phagocytosis and cell 
trafficking. Besides actins and heat-shock proteins, the absence of 
key immune-regulatory molecules amongst the most-significantly 
up-regulated genes in IB remains to be investigated; one possible 
explanation could be that, similarly to schistosome parasites 



infecting susceptible Biomphalaria spp. [14], 0. viverrini manipulates 
the innate defense system of susceptible Bithynia, effectively 
resulting into the inability of the hemocytes to recognize the 
parasite and/or into the suppression of the snail humoral 
responses against the parasite invasion. This hypothesis is also 
supported by the fact that, out of the DEGs identified in this study, 
the vast majority (i.e. ~90%) were molecules whose levels of 
transcription were up-regulated in UB, thus reinforcing a potential 
role of 0. viverrini in mechanisms of immune suppression in 
susceptible snails [14]. 

Taken together, the similarities in gene expression profiling 
between parasite-free specimens of wild-caught B. siamensis 
goniomphalos examined in the present study and schistosome- 
resistant populations of B. glabrata, and between infected Bithynia 
and schistosome-susceptible Biomphalaria, lead to the speculation 
that an O. viverrini-resistant population of Bithynia snails may occur 
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Figure 2. (a) Enriched Gene Ontology (GO) terms assigned to differentially expressed genes (DEGs) in uninfected and Opisthorchis 
wVe/r/nAinfected Bithynia siamensis goniomphalos snails (UB and IB, respectively), according to the categories 'biological process', 
'cellular component' and molecular function; (b) enriched KEGG biological pathways assigned to DEGs in UB and IB, respectively. 

doi:1 0.1 371 /journal.pntd.0002765.g002 



in natural freshwater bodies of Northeast Thailand. However, it 
must be pointed out that the snail specimens herein examined 
were screened for the presence of 0. viverrini DNA only; therefore, 
the occurrence of underlying infections by other snail pathogens in 
both UB and IB cannot be excluded. Ultimately, in order to test 
the hypothesis that O. viverrini-resislaxA B. siamensis goniomphalos may 
occur in Thailand, the development of adequate and efficient 
methods to successfully rear and maintain populations of this snail 
species under laboratory conditions is necessary. Once an optimal 
rearing protocol is established, future studies could investigate 
changes in gene expression of key genes at different time points 
following experimental infection of Bithynia specimens with 0. 
viverrini; such studies will allow to further explore the role/ s of 
selected molecules, such as proteases and actins, in the cascades of 
molecular events that determine susceptibility or resistance of 
snails to trematode parasites. 



A large collaborative project aimed at generating the first 
reference whole genome sequence of B. glabrata is currently 
underway (http://biology.unm.edu/biomphalaria-genome/). Once 
completed, the B. glabrata genome will represent an unprecedented 
resource for in-depth studies of snail- trematode interactions, as well 
as for comparative analyses between and among species acting as 
intermediate hosts of parasitic flukes, and for the discovery of the 
molecular and evolutionary mechanisms associated to the selection 
of parasite-resistant snail populations [67,68] . Similarly, the wealth 
of molecular data that is currently available for 0. viverrini 
[10,1 1,69], together with future efforts to define its entire genome 
sequence, will set the basis for the identification of genes and gene 
products essential for the recognition and infection of susceptible 
snail hosts. Taken together, these advances will provide the scientific 
community with new tools for the development of strategies to 
control the spread of opisthorchiasis and Opistorchis induced-CCA. 
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Figure 3. The top ten up-regulated genes linked to biological pathways associated to immune responses [based on orthology to 
molecules in the Kyoto Encyclopedia of Genes and Genomes (KEGG)] in each uninfected and Opistorchis wVesr/nAinfected Bithynia 
siamensis goniomphalos snails (UB and IB, respectively). 

doi:10.1371/journal.pntd.0002765.g003 
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